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I. INTRODUCTION 


A. MOTIVATION 

Accurate modeling of progressive flooding is essential in assessing the survivability 
features of a given ship. A realistic flooding and counter flooding scenario is a very 
complex dynamic system including complicated vessel motions, time-varying external 
excitation, changing hull mass and hydrodynamic properties, progressive flooding 
through damage openings, internal sloshing of flooded water, and de-watering control 
actions. Although, due to the complexity of the problem, the traditional method of 
approach 1s through large amplitude static analyses, the need for a more accurate analysis 
has been recognized. In order to bridge this gap, recently there has been a large amount of 
research focusing on dynamic stability prediction taking into consideration water-on-deck 
and large amplitude time-dependent motions and excitation [Ref. 1]. Equally important, 
however, is the development of a design methodology for incorporating elements of 
man/machine interface into the problem; in other words effects of decision-making and 
control strategies. This presents an additional level of difficulty since combinations of 
discrete events and continuous differential equations describe such systems. The long 
term goal that is initiated with this thesis is to study the development of an improved 
damage stability analysis and design tool, which will address the effects of control 


decisions and vessel dynamics in a progressive flooding scenario. 


B. BACKGROUND 

When a ship suffers underwater damage, there is a rapid influx of flooding water, 
usually followed by a period of slower progressive flooding. Progressive flooding is 
defined as flooding which extends beyond the compartments that are open directly to the 
sea. This results in flooding of compartments whose hull boundaries, but not interior 
bulkheads, are still intact. Existing computer models at the Naval Postgraduate School 
utilize the FORTRAN computer language and formatted input files to model progressive 
flooding [Ref. 2]. Any change in model design conditions, flooding damage, or damage 
control abilities requires a completely new FORTRAN input file. The purpose of this 
thesis is to include the capabilities of the original computer progressive flooding model in 
a more modular and expandable format, using the user-friendly MATLAB and 


SIMULINK programs. 


C. DEFINING THE MODEL 

1. Hole Model 

Holes can be modeled as a short-tube orifice as shown in figure | with a diameter 
(hole size) much larger than its tube length (ship’s skin thickness). The flow discharge 
coefficient for this case is 0.816 [Ref.3: pp.42-43]. 

2. Hull Model 

The required hull characteristics are determined using the same table of offsets 
used by Dawson (Hull.dat) [Ref. 2: pp. 127-143] and MATLAB programs arvol .m and 
drafts.m tocalculate the locations of center of buoyancy and center of floatation, 


waterplane area, and waterplane moments used in the program equations. Both Matlab 


programs are included in the appendix. For the FFG-7 hull model used here the 
following values were calculated: 
Table 1: Hull Model Characteristics 


Center of Buoyancy 


>.C.61 | Caablep least ee ea 209.9 ft aft of FP 
CE) CEs ee eS AR ea 0.0 ft from CL 
VAC) CRM ce tr eee 10.52 ft from baseline 


Nelo eo oc ee 227i att Ol FP 

XD al. . Ae eae. 0.0 ft from CL 
Wyiaterplane arCammmmensreere meee ccc, 14030 ft * 
Longitudinal Moment of Inertia............ Posse tt” 
Transverse Moment of Inertia............... 1.8148 x10° ft* 


3. Equilibrium Model 

The added weight method [Ref.4: pp.76-79] is used within the program to 
determine the equilibrium condition (displacement, heel, and trim) of the hull model for 
any state of flooding. This method assumes hull characteristics are constant and is 
accurate for “small angles”. The added weight method treats flooding water as added 
weight and treats the hull as if it were intact. The effect of the added weight is divided 
into a weight, which causes the ship to immerse evenly (parallel sinkage) and moments 
which cause the ship to list and trim about its center of floatation. The small angle 
assumption assumes that the ship lists and trims about a fixed center of floatation. In 


reality, since the hull is not cylindrical, the location of the center of floatation will change 


with changing list and trim. This method allows much simplified calculation of hull 
characteristics compared with the lost buoyancy method. Since the weight added is liquid 
water which is free to move as the ship moves (and thus does not contribute to pitch and 
roll resistance), a free surface correction is added to the longitudinal and transverse 
moments of inertia to account for the reduced resistance of the hull to pitch and roll. 

4. Compartment Model 

The ship compartments are modeled as rectangular boxes. Extension to more 
realistic shapes is one of the recommendations for further expansions to the capabilities 


of the program. 


D. SIMULINK 

The SIMULINK program provides a graphical representation of mathematical 
relationships and equations to enable the user to determine the response of system 
odels. Graphical symbols and connections are used instead of numbers and 
mathematical symbols, allowing more intuitive use and understanding of complex 
mathematical relationships between system components. The SIMULINK program also 
allows easy modification of input parameters and import and export of data from 
MATLAB programs. Lastly, SIMULINK allows the user to observe system or individual 


component behavior during the simulation through the use of output graphs. 


E. OBJECTIVES 
The objective of this thesis is to develop a progressive flooding model to utilize the 


superior user interface and data exchange capabilities of the SIMULINK program. The 


program must be able to accurately calculate flooding rates and the subsequent response 


of the ship. The program must also be user-friendly, easily modified and have the 


potential for further expansion by subsequent users. To achieve these goals the following 


programs were written by the author: 


2. 


1. Programs 


a. Pflood.m 

Pflood.m is the main program created by the author. This program calculates the 
flooding and ship response caused by a hole of user-chosen size and location in the 
ship’s hull. The flooding rate, flooded volume, heel and trim are calculated and 
updated each second. In addition, the user can specify the number, time of 
activation and size of pumps used for dewatering and the duration of their use. 

b. Arvol.m 

Arvol.m is a MATLAB program written by the author to determine the center of 
buoyancy of the hull model defined by the Hull.dat data file using Simpson’s rule 
integration techniques. 

c. Drafts.m 

Drafts.m is a MATLAB program written by the author to determine the waterplane 
area, center of floatation, and second moment of inertia of the hull waterplane for 
a given draft. 

Validation 


Program performance was validated by comparing the program output to expected 


ship behavior for a given flooding condition. 





Il. SIMULINK PROGRAM 
The schematic representation of Pflood.m is shown in Figure |. For illustrative purposes, 
Pflood.m is shown in its “expanded” form. By leaving Pflood.m in expanded form, any 
function or variable can be observed or changed by the user. For ordinary use, the various 
subroutines would be grouped into subsystem blocks with only the user-defined variables 
and the outputs represented by discrete function blocks. Placing the program in 
subsystem block form has the added advantage of greatly enhancing program readability 
and user-friendliness. Specific details of the individual block diagram elements are 


included in Chapter II. 
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1: Pflood.m Simulink Program 


Figure 


Iii. PROGRAM AND SUBROUTINES 


A. FLOW RATE THROUGH HOLKE(S) 

The flow through the damage, either primary or secondary, can be modeled as 
turbulent flow through an orifice of negligible length, or a hole. Flow rate through a hole 
is a function of the hole size (cross sectional area), hole shape, and pressure difference 
across the hole. The hole size and shape are related to the discharge coefficient, Cg. A Cg 
of 0.816 is used in this application in order to remain consistent with Dawson’s previous 
computer flooding model. This number is also consistent with theoretical predictions of 
flow through orifices as explained in the previous chapter [Ref. 3: pp. 42-43]. The 
pressure difference is determined by the difference in height of the liquid on either side of 


the hole. The equation is 


O=CdxAx.J2g x(h2—-h) 


where Q is the volumetric flow rate, A is the cross sectional area, Cd is the discharge 
coefficient, g is the acceleration due to gravity, and ho-h; is the liquid height difference 
across the hole. This height difference is a function of the level of the flooding water 
inside the flooded tanks, as well as the draft and trim of the ship. The latter quantities 
will have to be calculated using the added weight method as explained in the previous 
chapter. The Simulink representation of the above equation is shown in Figure 2. 

The program allows continuous adjustment of the values of hl and h2 to account for 


increased immersion of the exterior hole and increased compartment water height. 





Height of Wate 


gure 2: SIMULINK Hole Flow Symbology 
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B. NET FLOW INTO A COMPARTMENT 

Flooding water removed from the compartment by pump activation or by flow 
into adjacent compartments is subtracted from the water flow rate into the compartment. 
Two way flow due to water removal from one of two connected flooded compartments is 
also possible, and is accounted for by taking the square root of absolute value of the 
liquid height difference and then multiplying the result by the sign of the height 
difference. This method avoids program failure due to trying to take the square root of a 
negative number. This program also continuously updates the values of h1 and h2 to 
account for increased immersion of the exterior hole and increased compartment water 


height. 


C. HEIGHT OF WATER (HOW) 

Once the net volumetric flow rate of water is known, integrating the flow rate 
results in the net volume of water added to the compartment. This volume, divided by the 
horizontal area of the compartment, and divided again by the permeability of the 
compartment, gives the height of water in the compartment. A saturation function ensures 
that the height of water in the compartment is always between zero and the maximum 


height of the compartment. Figure 3 shows the Simulink height of water function. 


D. MOMENTS DUE TO ADDED WATER 


The added weight method treats the added weight as a combination of a weight 


which causes the ship to sink evenly (parallel sinkage) and transverse and longitudinal 
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Figure 3: SIMULINK Height of Water Terminology 


moments which cause the ship to heel and trim respectively. Multiplying the height of 
water (HOW) by the product of the weight density of the liquid ( pg ), the compartment 
permeability (per), and the compartment horizontal area (A) determines the weight added 
to the compartment. Multiplying the added weight by the moment arm (Zirrans/iong) from 
the ship’s center of floatation to the center of gravity of the added weight gives the 
moment about the center of floatation in the longitudinal and transverse directions. In 


equation form: 
Moment = (pgAX HOW X per) X Zarans / tong 


The Simulink representation of the above equation is shown in Figure 4. 


E. KG (VERTICAL POSITION OF CENTER OF GRAVITY) 

KG (Vertical Position of Center of Gravity) is determined by multiplying the 
original weight of the ship by the original KG, adding the product of the weight of the 
added water and the height of the added water center of gravity above the keel, and 


dividing the sum by the new ship weight. The equation for KG is: 


KGnew = (KGotp X Wot + KGapp X Wapp) + Wnew 


The Simulink representation is shown 1n Figure 5. 


F. VERTICAL CENTER OF BUOYANCY (KB) 
The ship’s new vertical center of buoyancy (KB) is determined by multiplying the 
original KB by the original displacement, adding the product of the added weight and the 


original draft plus one half of the parallel sinkage, then dividing the product by the new 
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Figure 4: SIMULINK Moment Terminology 
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Figure 5: SIMULINK KG Terminology 


weight. In equation form: 

K Bnew = ( KBoric X Woric + KBapp X Wapp) + Wnew 
This 1s based on the assumption that the hull 1s effectively wall-sided, and that the vertical 
center of buoyancy of the newly submerged portion is located at one half the added draft. 
Both of these assumptions are consistent with the small angle assumption that was 
mentioned in the previous chapter. If the small angle assumption is no longer true, then 
the true position of the center of buoyancy will have to be calculated. This can be done by 
employing numerical integrations and utilizing the provided table of offsets, or by 
interfacing the program with a standard external ship hydrostatics calculations program, 
such as SHCP [Ref. 2: pp. 24-25]. Such an interface is relatively straight forward since 
the computational environment of Matlab/Simulink already allows incorporation of 


external ASCII files. Figure 6 shows the Simulink representation of the KB equation. 


G. HEIGHT OF THE METACENTER ABOVE THE CENTER OF BUOYANCY 
(BM) 


The height of the metacenter above the center of buoyancy (BM) 1s determined by 
the second moment of the waterplane (about the waterplane centerline or longitudinal 
center of floatation) divided by the hull displacement. For the added weight method, a 
free surface correction is made by subtracting the second moment of the flooded 
compartment area from the waterplane moment. Typically, free surface corrections are 
done for liquids when calculating the final metacentric height. In such a case, the free 
surface correction appears as a reduction in the metacentric height or, equivalently, an 


increase in the virtual center of gravity. The same effect can be observed by modifying 
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Figure 6 : SIMULINK KB Terminology 


the second moment of the waterplane area in the calculation of the metacentric radius, 
BM. This is because the added weight in this case 1s seawater. With this modification in 
BM, the added weight method produces identical results with the lost buoyancy 
calculations (reference a naval architecture book). The BM Simulink terminology is 


shown in Figure 7. 


H. RIGHTING ARM (GM) 
The righting arm (GM) is determined by GM = KB + BM - KG. This equation is 


graphically represented by the triple summing point shown in Figure 8. 


I. HEEL AND TRIM ANGLE 

The angle of heel is determined by 

sin g = Volume X Zrrans + (Displacement x GM ) 

as shown in Figure 9. The total trim is determined by dividing the total trimming moment 
by the moment to change trim (MCT). MCT is equal to the total weight times the 
longitudinal GM, divided by the ship length at the waterline. The tangent of the trim 
angle is equal to the total trim divided by the ship’s length at the waterline. If the large 
angle assumption needs to be utilized, the moment to change trim will have to be 
imported as an external array. This is not expected to cause significant difficulties since 
MCT is readily available from the hydrostatic diagrams of a given ship. The Simulink 


trim angle terminology is shown in Figure 10. 
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Figure 7: SIMULINK BM Terminology 
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Figure 8: SIMULINK GM Terminology 
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Figure 9: SIMULINK Heel Angle Terminology 
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Figure 10: SIMULINK Trim Angle Terminology 


J. ADDITIONAL SUBMERGENCE 

The additional submergence of the hole due to list 1s determined by multiplying the 
sine of the heel angle by the transverse distance from the transverse center of floatation to 
the hole. Similarly, the increased hole depth due to trim is obtained by multiplying the 
sine of the trim angle by the longitudinal distance to the center of floatation. The new 
hole depth is the sum of the original hole depth and the depth increases due to heel and 
trim. In algebraic form, these relations are shown by 


h2new = h2oww + (ZTRANS , hole X SiN @) + (ZLONG , hole X Sin( trim )) 


and in graphical form they are represented by the diagram of Figure |]. 


K. PROGRAM USAGE 
To use Pflood.m, the user needs only select “simulation” from the SIMULINK 
menu, then select “start”. Modifications to the simulation parameters or the program are 
~ easily accomplished by “double-clicking” on the symbol of interest. Additionally, 
SIMULINK allows the user to easily add output displays in order to observe the behavior 


of any desired variable during the simulation. 
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Figure 11: SIMULINK Hole Submergence Terminology 


IV. RESULTS 


Results for a test case consisting of a compartment with a 4 square foot hull 
penetration (compartment 1) and a 4 square inch watertight bulkhead penetration 
(compartment 2). Pumps were activated at 10 minutes (600 seconds) which completely 
dewatered compartment | in order to test the ability of the program to account for 
backflow through holes. As shown in figure 12, the four square foot hull penetration 
completely floods compartment | essentially instantly. Activating the pumps 
dewaters compartment 1, instantly removing the water entering through the hull 
penetration. Similar results could have been achieved by changing the characteristics of 
the hull penetration to simulate hull repairs. 

Figure 13 illustrates how the program accurately predicts backflow into 
compartment! from compartment 2, which had no pumps assigned. A saturation function 
ensures that backflow from compartment 2 will stop when compartment 2 water height 
drops to the level of the bulkhead penetration. 

The predicted trim response of the ship is shown in Figure14. The ship initially 
trims as a result of the rapid flooding of compartment |, then gradually increases its trim 
angle as compartment 2 floods. The instant dewatering of compartment | then causes the 
trim angle to decrease, followed by further trim angle reduction as compartment 2 is 
dewatered. Without altering the pump or hole configuration, the ship will assume a steady 
trim angle corresponding to compartment 2 flooded up to the height of the bulkhead 


penetration. 
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Figure 12: Height of Water in Compartment | 
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Figure 13: Height of Water in Compartment 2 
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Figure 14: Ship Trim Response 
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V. DISCUSSION, CONCLUSIONS, AND RECOMMENDATIONS 


A. DISCUSSION 

This thesis has successfully developed the foundation for a modular and 
expandable progressive flooding computer model. The computer model accurately 
predicts the quasi-static ship response to progressive flooding and is easily modified by 


the user to explore ship response to any desired damage situations. 


B. CONCLUSIONS 
Based on the above results, the following conclusions are drawn: 


1. Pflood.m is an effective design tool for modeling and studying the effects of 
progressive flooding. 


2.  Pflood.m is easily modified to model virtually any damage scenario. 


3. Pflood.m can be expanded to include different ship progressive flooding and 
equilibrium models. 


C. RECOMMENDATIONS 
The following items are recommended to further expand Pflood.m and enhance 
its utility as a design tool: 
Incorporate a hull characteristics look-up table into the Pflood.m program 
to expand the ship response modeling ability beyond “small angles” of 


heel and trim and non-rectangular compartments. 


2. Explore the effect of pump configuration on dewatering effectiveness and 
ship response. 
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APPENDIX B. ARVOL.M 


arvol.m 

by Russell wright 

calculates station areas and center 0’ bouvancy 
EOmn =16.9 £6 


dP dP dP a? 


lead hulll .dat 
=zeros (41,3); 
Zaza 16.9 

Bon y=1l:41 


beg=23* (1-1) +1; 
endd=23 <i. 
a(i,1)=hulll (beg) ; 
z=nulli (beg: endd,3): 
eee o> Zz 
ae ee ee Oa a) = 

else 

%* y=hulll (beg:endd,2); 


zp=hulll (beg:endd,3); 
Wawel (6. 9=2oqt jor 22:16.9; 
yvy=hulll (beg: endd, 2) ; 
y=spline(zp,yy,2Z); 
mezi2)-z Ly; 

hhh (i) =h; 
mour = /020-rr=0: 
two=0.0;to=0; 


for j= L:121 
EOUT=fOur+y\ ce 30 
Pretaec(2 3) Vi2en. 

end 


zor k=1-10 
two= twory(2*k+1); 
PO=bOry( 2 Kp) Zz (27k+1); 
end 


% calculate station area using simpson’s rule: 
eu 2) =2~(y(1)+4* tour+2 *twoty(23)) *n/3- 

% calculate station centroid (z bar) using Simpson’s rule: 
eee =o iy (1) 2 (1) +4 eee? * cory (23) *2(23)) *h/(3"ali,2)): 


end 
end 


aa=zeros(101,3); 
bnm=((a(37.1)-a(3,1))/100); 

peewee, |) =(a(3,1):nh:a(37,1))‘; 

eeae., 2) =spline(a(3:37,1),a(3:37.2),.aa(:,1)); 
emer s | —SD line (a(3:37,1),a(2-37,3),aa(:.1)); 


meme=0;too=0;7 fre=0; tw=0; ro=0; tu=0; 
fom 2=1:50 
fore=fore+aa(2*1,2); 
mee=Ere+aa(2*i,2)*aa(2*1,3); 
moa—Lotaa(2*1,2) *aa(2*1,1); 

end 


moe i=1:49 
too=too+aa(2*i+l,2); 
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Beet oa er ely, 2)aa (2 <2 <1, tie 
Beeettea( evi), 2) *aa(o*i+i. ie 
end 


vol=(hh/3)*(aa(1,2)+4*fore+2*tootaa(101,2))- 


Spee a eT ee age NPG Ay em lay) se ooo, 
eee ey) 4 fond stu aa (101 rayeeenaion )) neene, 
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APPENDIX C. DRAFTS.M 


oP 
foo) 
dP 
dP 
ae 
oe 
dP 
of 
ae 
oe 
de 
dP 
of 
ae 
dP 
oP? 
of 
of 
uf 
cP 
ae 
dP 
dP 


SSSSESEESELSESSESTESSSSFESELEEEEEEESESEEESEEEEESYSEE SESS 


of of 


erate sZom 
BY xzsussell wright 
calculates wp points for a given draft 


oP oP oP of of cP 


load hullidat 


ind=0; 
Vincl=zeros (417377, 
22= Loe; 


Zine 2s 

ZEneZ=2z 

Pemer =. +41 

a 
beg=23* (1-1) +1; 
ence=25 "2 
wpola(i,1)=hulll (beg, 1) ; 
Z=ni iT ibeg2enda:.s )-; 
V=ebulii(Seq:enda, 2); 


Wola te) =interol(2,V,22, seline ); 
ae aes) oy Ee 
We lanCi.2) 
yinel(i, 
Vale is. 
wola(l, i 
end 


i—— i 


Vine Were) =e seq sl ).- 

ViEnet 2 Inte sol 2, Veacine |) Spline’); 
Pe Vincl( iy) 2 jeeco 0 

ViEtiC Lilie =o aU 

end 

Vine l( i338 )=22; 


1f£ wpla(i,2)==0 


)= 
Vite dt. 2) Go; 
Yarnell 3) =0; 
end 
end 


Vine l=V ine s87 2): 


KEV Ine 

Mane io) eee Sl poeme (1) )/34evinel (35,1). 
VV =eVINC L222; 

y=spline(xx,yy,x); 


h=x(2)-x(1); 
Jedelotl(eo) = ole 
fou = 702 0- EH=0.- 


two=0.0;to=0; 
for ==) ely 


POUT = cour — (2-4). 
Eea=Er res a) 02" 4) 5 
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end 


Or t= the 
two= twory(2*k+1); 
Ceo=corvy 2*k+l)=x(2 k+l) 
end 


% calculate wp 1/2 area using simpson’s rule: 
aZg=2 "(vi 1) +4* Sour-2*two+ry( 35) )*h/3- 

% calculate wo centroid (x bar) using simpson’s rule: 
MDar=2 ty (1) 50) +4“ tr+2* tory (35) * oo ny C3 al) 


de 


calculate long. moment oo’ inertia about xbar 


Gee) =) den 


Prateek (2g) yo" 5) 
end 


for k=l=]6 


Eeo=stoty( 27 k+1)*(x(2*k+1) 72): 
end 


= 2 sy Cec) ta * Ere * Coty (35) > €x(35) > 2). by 3: 


Lsubl=1o0-a2*xbar%*2 
% calculate trans moment o’ inertia 


ssp Be 
Eo=0- 
CoE wa lols, 


rere ty 2 9) 3): 
end 


for k=1:16 


EO=bOf (712 kr) 3) 
end 


lsubt=2* (ly (1) °3) et er+2*to+(y (35) 73) ) ag, 
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